Stable liquid Hydrogen at high pressure by a novel ab-initio molecular dynamics 
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We introduce an efficient scheme for the molecular dynamics of electronic systems by means of 
quantum Monte Carlo. The evaluation of the (Born-Oppenheimer) forces acting on the ionic posi- 
tions is achieved by two main ingredients: i) the forces are computed with finite and small variance, 
which allows the simulation of a large number of atoms, ii) the statistical noise corresponding to the 
forces is used to drive the dynamics at finite temperature by means of an appropriate Langevin dy- 
namics. A first application to the high-density phase of Hydrogen is given, supporting the stability 
of the liquid phase at ~ SOOGPa and ~ 400iC. 

PACS numbers: 47.11. Mn, 02.70.Ss, 61.20.Ja, 62.50.-|-p 



The phase diagram of Hydrogen at high pressure is still 
under intense study from the experimental and theoret- 
ical point of view. In particular in the low temperature 
high-pressure regime there is yet no clear evidence of a 
metallic atomic solid, and either a molecular solid phase 
can be favoured even in this regime^, or the liquid phase 
can be stabilized at low temperature by increasing the 
pressure (see Ref.([l[)). 

Indeed, for high pressures around SOOGFa, a two 
fluid (proton and electron) superconducting phase, in- 
duced by the strong electron-phonon coupling, has been 
conjectured and unusual quantum properties have 
been later predicted^. In this work we use an improved 
afe-initio molecular dynamics (AMD) by using accurate 
forces computed by Quantum Monte Carlo (QMC). We 
present preliminary results, showing that the liquid phase 
is energetically stable, due to the strong electron correla- 
tion, at least within the Resonating Valence Bond (RVB) 
variational approach [Bj, which is very accurate also in the 
solid phase. 

AMD is well established as a powerful tool to investi- 
gate many-body condensed matter systems. Indeed, pre- 
vious attempts to apply Quantum Monte Carlo (QMC) 
for the dynamics of ions 6] or for their thermodynamic 
properties are known, but they were limited to small 
number N of electrons or to total energy corrections of 
the AMD trajectories, namely without the explicit calcu- 
lation of the forces. Indeed, the technical achievements 
that we are going to present in this letter are particu- 
larly important for the simulation of liquid or disordered 
phases by QMC. 

Calculation of forces with finite variance. The simplest 
method for accurate calculations within QMC, is given by 
the so called variational Monte Carlo (VMC), which al- 
lows to compute the variational energy expectation value 
EvMC = of a highly accurate correlated wave 

function (WF) by means of a statistical approach: 
electronic configurations {a;}, with given electron posi- 
tions Ti and spins ct; — ±1/2 for i = 1, • • • A^, are usually 



generated by the Metropolis algorithm according to the 
probability density oc iprix)'^- Then Evmc is com- 
puted by averaging statistically over ji^ the so called local 
energy etix) = ^tli^Jf , namely Evmc = J dfix eLix), 



where / d/ix indicates conventionally the 3A^ multidimen- 
sional integral over the electronic coordinates weighted 
by tp^{x). In the present work we assume that the WF 
iPt{x) = {x\iPt) = J X detA is given by a correlated 
Jastrow factor J times a determinant D oi a. N x N ma- 
trix A, such as for instance a Slater determinant. The 
main ideas of this approach can be straightforwardly gen- 
eralized to more complicated and more accurate WF's, 
as well as to projection QMC methodsjsl more accurate 
than VMC. 

The eflticient calculation of the energy derivatives. 



namely the forces /s 
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■Na, 



'R. dR~ 
where Na is the number of atoms, is the most impor- 
tant ingredient for the AMD. Within VMC they can be 
computed by simple differentiation of EvmCi using that 
not only the Hamiltonian H but also ipT depend explic- 
itly on the atomic positions Ri. This leads to two dif- 

/- , the 

■' Ri' 



ferent contributions to the force fs — f - 

•' Ri ■' R 

Hellmann-Feynman / 



HF 
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and the Pulay one , where: 



(1) 



d^x {eL{x) - EvMc)dj} log\il;T{x) \ (2) 



However in order to obtain a statistically meaningful 
average, namely with finite variance, some manipula- 
tions are necessary because the first integrand may di- 
verge when the minimum electron-atom distance van- 
ishes, whereas the second integrand is analogously un- 
bounded when an electronic configuration x approaches 
the nodal surface determined by ipTix) — 0. By defining 
with d (6) the distance of x from the nodal region (the 
minimum electron-atom distance), CLix), dj^ logiprix) — 
1/d ( {x\d^,H\x) ~ 1/(5^), whereas fix — d"^ {fix — ^^), 
leading to an unbounded integral of the square integrand 
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in Eq.([2l) (EqlT]), namely to infinite variance. The in- 
finite variance problem in Eq.([T]) was solved in several 
ways. Here we adopt a very elegant and efficient scheme 
proposed by CafFarel and AssarafQ. Instead the infinite 
variance problem in Eq.([2]) was not considered so far, 
and this is clearly a problem for a meaningful definition 
of ionic AMD consistent with QMC forces. 

In this letter we solve this problem in the following 
simple way, by using the so called re-weighting method. 
We use a different probability distribution /j,^ c>c Tpci^)'^! 
determined by a guiding function ^g{x): 



(3) 



where R{x) oc iprix) ^ for d — > is a "measure" 
of the distance from the nodal surface iPt{x) — 0. By 
assumption ipx may vanish only when det A — (J > 0) 
and therefore R{x) is chosen to depend only on A. For 
reasons that will become clear later on we have adopted 
the following expression: 



R{x) = 1/, 



N 



(4) 



Then the guiding function is defined by properly regular- 
izing R{x), namely: 



R'ix) 



R{x) if R{x) > e 

e(i?(x)/e)^(^)A if rIx) < e 



(5) 



The non obvious regularization for R{x) < e instead of 
e.g. R'^{x) = Max[e, R{x)] was considered in order to sat- 
isfy the continuity of the first derivative of ijjGix) when 
R{x) = e, thus ensuring that '4'g{x) remains as close 
as possible to the trial function Vt- In this way the 
Metropolis algorithm can be applied for generating con- 
figurations according to a slightly different probability 
H'^ix) and the exact expression of can be obtained by 
the so called umbrella average: 

-2jdnl S{x){eL{x) - EvMc)djiJogijjT{x) 



Now, the re- weighting factor S{x) = {^t{x) / ipaix))^ — 
Min [1, (i?(2;)/e)]^~^^*-'^-'^^ oc d? , cancels out the diver- 
gence of the integrand, that was instead present in Eq.(l2]). 
Hence the mentioned integrands in the numerator and 
S{x) ( < 1) in the denominator of Eq.(l6|) represent 
bounded random variables and have obviously finite vari- 
ance. In this way the problem of infinite variance is def- 
initely solved within this simple re-weighting scheme. 

We show in Fig.((T]) the efhciency of the method for 
computing the Pulay force component acting on a Hy- 
drogen proton at = 1.31 in a bcc lattice. As it is clear 
in the plot for the N = 128 case, the difference between 
a method with finite variance and the standard one with 



infinite variance is evident. In this way the 3iV"^ x ?>Na 
correlation matrix aqMC, defining the statistical corre- 
lation between the force components, can be efficiently 
evaluated: 

aQMc{R) -< (4- < 4 >)^h- < 4 >) > 

where the brackets <> indicate the statistical average 
over the QMC samples. The correlation matrix aqMCi 
that within the conventional method is not even defined, 
will be a fundamental ingredient for a consistent AMD 
with QMC forces and therefore the solution of the infinite 
variance problem is crucial for this purpose. 
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FIG. 1: (color online). Evolution of the integrand in Eq.(|6]) 
as a function of the Monte Carlo iterations. Each new sample 
is obtained after 2A'' Metropolis trials. 



Langevin dynamics. In the following derivation we as- 
sume that ions have unit mass, that can be generally ob- 
tained by e.g. a simple rescaling of lengths for each ion 
independently. For clarity and compactness of notations, 
we also omit the ionic subindices i when not explicitly 
necessary. Moreover matrices (vectors) are indicated by 
a bar (arrow) over the corresponding symbols, and the 
matrix-vector product is also implicitly understood. We 
start therefore by the following AMD equations for the 
ion coordinates R and velocities v: 



V = -jiR)v + fiR)+rf{t) 
R = V 



(8) 
(9) 



By using the fluctuation-dissipation theorem the friction 
matrix 7 is related to the temperature T (henceforth the 
Boltzmann constant ks — 1) by: 



7(i?) = ^aiR) 



(10) 



where a{R) is generally a symmetric correlation matrix: 
<ri,it)rf,{t')>^6it-t')a{R). (11) 
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It is important to emphasize that, as a remarkable gener- 
alization of the standard AMD used in [lO|] , in the present 
approach the friction matrix 7, may depend explicitly on 
the ion positions R, so that Eq. ([TU]) can be satisfied even 
for a generic correlation matrix a{R). In fact we have 
the freedom to consider a QMC contribution in a{R): 



a{R) = ao + Aq aQMciR) 



(12) 



where Aq > and ao is the identity matrix / up to 
another positive constant ao, ao = ctol, and aQMc{R) 
can be estimated by Eq. ([7]) • In the following we will show 
that, for appropriate ao, Ao > 0, it is possible to follow 
the Langevin dynamics by means of noisy QMC forces. 

Integration of the Langevin dynamics. Henceforth the 
velocities Vn are computed at half- integer times i„ — 
whereas coordinates i?„ are assumed to be defined at 
integer times Rn = R{tn)- Then, in the interval tn — 
^ < t < tn + ^ and for At small, the positions R are 
changing a little and, within a good approximation, the R 
dependence in the Eq.([S]) can be neglected, so that this 
differential equation becomes linear and can be solved 
explicitly. The closed solution can be recasted in the 
following useful form, where the force components appear 
corrected by appropriate noisy vectors t): 



Vn+l 
Rn+1 
f 



e-^^*^/„ + f(/(i?„)+f)) 



(13) 
(14) 
(15) 



V = 



Rn + AtVn+l + 0{At^) 

^ J dte^'^'-'"^rf{t) (16) 



2sinh(4^7) 



By using that a — 2Tj from Eq. (fTO|) . and that its depen- 
dence on R can be consistently neglected in this small 
time interval, the correlator defining the discrete (time 
integrated) noise ff can be computed explicitly: 



sinh(Ai7) 
4sinh(^7)2 



a 



(17) 



This means that the QMC noise has to be corrected in a 
non trivial way as explained in the following. 

Noise correction. The QMC noise is given during the 
simulation, and therefore in order to follow the correct 
dynamics another noise 77 '^^* has to be added to the noisy 
force components in a way that the total integrated noise 
is the correct expression (ITTl) . i.e. i] = ff'^^'^ + f]QMC- By 
using that the QMC noise in Eq.© is obviously inde- 
pendent of the external noise, we easily obtain the corre- 
sponding correlation matrix: 



TABLE I: Comparison of the total energy per proton 
(Hartree) for Hydrogen in the bcc lattice at Vs = 1.31 com- 
pared with the published ones with lowest energy (to our 
knowledge) . 



N 


EvMC /Na 


Evmc/Na\U] 


Edmc/Na 


Edmc/Na\U] 


16 


-0.48875(5) 


-0.4878(1) 


-0.49164(4) 


-0.4905(1) 


54 


-0.53573(2) 


-0.5353(2) 


-0.53805(4) 


-0.5390(5) 


128 


-0.49495(1) 


-0.4947(2) 


-0.49661(3) 


-0.4978(4) 


250 


-0.49740(2) 




-0.49923(2) 




432 


-0.49943(3) 








00 


-0.501(1) 




-0.503(1) 





On the other hand, after substituting the expression ([12 
in Eq. pH)) 7 = 2^(0^0 + Ao aQMc) and using the ex- 
pression (jl7p for a', we obtain a positive definite matrix 



in Eq. (fT8|) for At < Aq 11|. Hence ff^^* is a generic 



a - OLQMC 



(18) 



Gaussian correlated noise that can be easily sampled by 
standard algorithms. After that the random vector 77'^^* 
is added to the force / + t]qmc obtained by QMC, and 
replaces / + 77 in Eq. (fT3l) . This finally allows to obtain 
an accurate AMD with a corresponding small time step 
error. The main advantage of this technique is that, at 
each iteration, by means of Eq. (fTU)) , the statistical noise 
on the total energy and forces (see Fig[21) can be much 
larger than the target temperature T, and this allows to 
improve dramatically the QMC efficiency. 

Optimization of the WF. In the following examples we 
consider a cubic box with periodic boundary conditions, 
and use a variational WF J x det A that is able to provide 
a very accurate description of the correlation energy, due 
to a particularly efficient choice of the determinant factor, 
that allows to describe the RVB correlations 12. [H. The 
WF contains several variational parameters, indicated by 
a vector (3, that have to be consistently optimized during 
the AMD. The Jastrow factor J used here depends both 
on the charge and spin densities, and it is expanded in a 
localized atomic basis. As it is shown in the table, the 
accuracy of our WF is remarkable. Indeed the small dif- 
ference between the so called DMC -providing the lowest 
possible variational energy within the same nodal surface 
of ^T- and the VMC energies clearly supports the accu- 
racy of our calculation. The size effects are very large for 
the metal and have been estimated {N — 00) following 
Ref.ie. 

In order to optimize the WF we use the recent method 
introduced in RefllSl. devised here in an appropriate way 
to optimize a large number of parameters during the 
AMD simulation, as described in Reffisl. This allows 
to remain efficiently within the Born-Oppenheimer en- 
ergy surface each time the ionic positions are changed 
according to Eq. fT^ . 

Application to high-pressure Hydrogen. We show in 
Fig-@ the evolution of the internal energy and corre- 
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-0.48^ 




500 1000 1500 500 1000 1500 2000 
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FIG. 2: (color online). Evolution of the internal energy and 
temperature vs the AMD with QMC forces. At = Ao = 
1.036/s, ao = 0.7fcsTa.M.. Bottom: points represent instan- 
taneous temperatures estimated by the average kinetic en- 
ergy, lines represent the target temperatures. They should 
coincide on average for At —> 0. The average energy, pressure 
and temperature in the last 0.5ps are —0.51319 ± 0.00003// 
(-0.5127±0.0001i/) 364K±5K (364±10ii:) and 335±2GPa 
(394 ± 5GPa) for TV = 128 (iV = 16), respectively. At each 
iteration the statistical noise on the total energy is ~ SOOOA'. 
The proton-proton g{r), averaged at the lowest temperature, 
is shown in the inset. 



spending temperature as a function of time with the 
proposed AMD with QMC forces, starting from the bcc 
Hydrogen sohd at = 1.31, considered henceforth. Al- 
though we have not studied the possible stabihty of all 
other solid phases yet, for N = 128 also the simple hexag- 
onal structure (SH) melts. This already provides a clear 
support to the liquid phase because these two atomic 
solids are the most stable ones, so far proposed at zero 
temperature. The proton-proton correlation function ob- 
tained starting from the two atomic solids is shown in the 
inset. 

We performed a finite size scaling analysis based on the 
comparison of our QMC results with the LDA ones[l6j. 
We found that LDA favours atomic solid phases respect 
to QMC. For instance the internal energy difference be- 
tween the liquid phase and the BCC solid one is only 
-0.0004 Ha in LDA while in QMC is -0.011 Ha. On the 
other hand the molecular solid fmhpc-cfljj) is also simi- 



larly preferred to the atomic solid by QMC, but appears 
to have much larger zero point energy corrections com- 
pared to the liquid. [13] Indeed we have studied quantum 
effects on protons by using the Wigner-Kirkwood expan- 
sion on the free energy and obtained 8.5{2)mH/ proton 
for the liquid, namely a much smaller correction than 
the SH [12.3(2)miJ/proton] and the molecular solid 
[13.1(2)m7J/proton] ones. Finally even the more accu- 
rate DMC does not affect the liquid stability, because it 
lowers the internal energy of all the phases studied by 
about the same amount (2 — Sm/J/proton). 

In conclusion we have shown that it is possible to 
make a realistic and accurate AMD simulation with QMC 
forces. We found that the bcc and SH sohd structures 
appear clearly unstable even at low temperatures where 
a molecular liquid has much lower internal energy, and 
is further stabilized by considering quantum effects on 
protons. This important finding highlights the present 
QMC technique as a possible and accurate alternative to 
study phase diagrams of materials and as a benchmark 
for other approximate methods. 
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